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ABSTRACT® 

The plane-stress equations of linear elasticity are used in conjunction with 
those of the boundary element method to develop a novel curved, quadratic 
boundary element applicable to structures composed of anisotropic materials in a 
state of plane stress or plane strain. The curved boundary element is developed to 
solve two-dimensional, elastostatic problems of arbitrary shape, connectivity, and 
material type. As a result of the anisotropy, complex variables are employed in the 
fundamental solution derivations for a concentrated unit-magnitude force in an 
infinite elastic anisotropic medium. Once known, the fundamental solutions are 
evaluated numerically by using the known displacement and traction boundary 
values in an integral formulation with Gaussian quadrature. All the integral 
equations of the boundary element method are evaluated using one of two methods: 
either regular Gaussian quadrature or a combination of regular and logarithmic 
Gaussian quadrature. The regular Gaussian quadrature is used to evaluate most of 
the integrals along the boundary, and the combined scheme is employed for integrals 
that are singular. Individual element contributions are assembled into the global 
matrices of the standard boundary element method, manipulated to form a system of 
linear equations, and the resulting system is solved. The interior displacements and 
stresses are found through a separate set of auxiliary equations that are derived using 
an Airy-type stress function in terms of complex variables. The capabilities and 
accuracy of this method are demonstrated for a laminated-composite plate with a 
central, elliptical cutout that is subjected to uniform tension along one of the straight 
edges of the plate. Comparison of the boundary element results for this problem 
with corresponding results from an analytical model show a difference of less than 
1 %. 
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INTRODUCTION 


Composite materials continue to see increased usage in space, aircraft, 
industrial, and recreational markets around the world. This upward trend is 
reflected in the recent advancements made for finite element methods; however, 
relatively few advances for boundary element methods (BEM) have been made for 
composite structures [1,2]. Boundary element advancements for composite 
structures are important, since boundary elements offer increased efficiency and 
accuracy over finite elements for infinite and semi-infinite regions, regions with 
large response gradients, and boundaries with complex geometry as well as reduced 
pre-processing time for discretizing a problem. 

The initial application of the BEM to bodies composed of anisotropic 
materials was made by Rizzo [3]. Cruse presented a further advancement of the 
BEM in 1971 that included a derivation of the traction and displacement 
fundamental solutions for an anisotropic plate as well as determining stress 
concentrations for anisotropic plates with circular and elliptical cutouts [4] . Recent 
texts written about the BEM have focused on basic problems associated with 
isotropic structures as well as advanced topics in the areas of fracture mechanics 
[5], acoustics, and the evaluation of complex integrals for body forces and three- 
dimensional volumes [6]. 

The primary contribution of the curved, quadratic boundary element 
developed in the present paper is to provide a capability for modeling anisotropic 
regions that have highly curved boundaries. Instead of modeling a circular or 
elliptical boundary by using many straight elements to approximate the geometry, a 
relatively few curved elements can be used to model the geometry exactly with no 
loss in solution accuracy. An important aspect of the curved, quadratic boundary 
element resulting from the material anisotropy is the treatment of complicated 
singular integrals. The integrand of the singular integrals is an elaborate expression 
in terms of complex variables that results from the anisotropic material behavior, 
quadratic shape functions, and Gaussian quadrature evaluation points. All integrals 
are evaluated using the quadratic shape functions and Gaussian quadrature, which 
provides very accurate modeling of curvilinear boundaries. Thus, by careful 
manipulation of the singular integrals a sophisticated scheme has been developed 
for analyzing structures that are composed of anisotropic materials, have complex 
curved boundaries, and are in a state of plane strain or plane stress. 

The objective of the present paper is to demonstrate a numerical method for 
evaluating elaborate, complex singular integrals that are used in the development of 
a curved, quadratic boundary element for structures composed of anisotropic 
materials that are in a state of plane stress or plane strain (referred to herein as plane 
anisotropic structures). This objective is accomplished by presenting the basic 
equations that govern a two-dimensional, linear elastostatic problem, the 
development of the boundary integral equation, the boundary element solution 
process, and the treatment of singular integrals using the developed method. A 
brief description of the modeling characteristics and analysis results for a classic 
structural mechanics problem is then presented to demonstrate the accuracy and 
computational efficiency of the curved boundary element for use with complex 
geometry. 



THEORY 


Boundary Integral Equation (BIE) Formulation 

First, the equations that govern two-dimensional linear elastostatics as given 
by Brebbia are stated [7]. The equations of equilibrium are written using indicial 
notation as, 


a kl i+bi=0 (k,l= 1,2) in Cl (1) 

where Ou are the components of the stress tensor, b; are the components of the 
body-force vector, and £2 is the problem domain enclosed by a boundary denoted by 
r. The equations derived in this paper follow the rules of indicial notation; 
specifically, partial differentiation is represented by terms containing a comma and 
repeated indices indicate terms that are summed. The kinematic relations are 

| (Uy + u j t i ) (i, j = 1 ,2) (2) 

where Ey are the components of the strain tensor and uj represents the components 
of the displacement vector. An additional equation that the strain components must 
satisfy when the displacement field is not the primary dependent variable is the 
compatibility equation. This equation in two dimensions is stated as 

E ll,22 +E 22,11 = ^ ^12,12 (3) 

Eq. 3 provides a necessary and sufficient condition for specified strain components 
to give displacements that are single-valued and continuous for simply connected 
regions. Since the problems of interest in the present paper are two-dimensional in 
nature, the constitutive equation for a state of plane stress is given for a plane 
anisotropic structure as 
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where Sq are the elements of the compliance matrix. In addition to eqs. 1-5, a set of 
boundary conditions is required to provide a complete description of the problem. 
The natural and essential boundary conditions are defined along the boundary T 
shown in Fig. 1 that consists of two distinct regions, Y 1 and IT. The tick marks 
along the boundary in Fig. 1 represent the separation between the two boundaries. 
On fi, displacements are prescribed by 


u;=Uj (on r 1) 


(6) 




Figure 1 : A general domain £2 with natural boundary conditions 
on IA and essential boundary conditions on Tj. 


and on Ta, the tractions are prescribed by 


Pk = a k/"i (onr 2 ) (7) 

where n z represents the components of the unit vector that is normal to T and that 

lies in the plane of the region £>. acted on by the traction. 

The boundary value problem is reduced to a single partial differential 
equation by expressing the stresses in eq. 4 in terms of the Airy stress function and 
then substituting the strains in eq. 4 into the compatibility equation, eq. 2. The 
general solution of the resulting partial differential equation is obtained by 
following Lekhnitskii’s method that uses a complex variable z given by z = x + py, 
where p is generally a complex number. Using this approach yields the following 
equation, known as the characteristic equation, which determines the values of p: 

^nP _ 7S 16 p +(2S 12 +S 66 )p — 2S 26 p + S 22 = 0 . (8) 

Lekhnitskii has shown that the four roots of eq. 8 are distinct and always imaginary 
as long as the material is not isotropic [8]. These roots can be denoted as, 

H= a k+^k P/fe = -ip* (At = 1,2) (9) 

where the lowercase letter i denotes V-T . Once the values for p are determined, the 
Airy stress function is expressed as a linear combination of complex potential 
functions. The forms of the complex potential functions are determined by the 
boundary conditions of a given problem. 

A wide variety of analytical and numerical techniques may be employed to 
solve problems of the type described above. The technique used here involves 
numerical integration of the boundary integral equation (BIE) that is derived below. 
Specifically, the BIE is used to relate known displacement and traction values on 
the boundary to unknown values of interest within the domain. A convenient point 
to begin the derivation of the BIE was developed by Cruse and others by using 



Somigliana’s identity which relates the value of the displacements at any point in Q. 
to the boundary values [9], u, and pj, such that 


u t (q) = J r u *j(%, x)pj (x)dr (x)-| r p*. (£, x)uj (x)dT( x) 


( 10 ) 


where p j is defined by the right hand side of eq. 7 and c, and x are functions of xi 
and X 2 . The values u* and p*- represent fundamental solutions that describe the 

displacements and tractions, respectively, in the j direction at point x corresponding 
to a unit-magnitude point force acting in the i direction of an infinite material 
applied at point The foundation for eq. 10 is the use of an arbitrary 
complimentary load case, represented by the terms with an asterisk. In this 
development, the load case is defined as a point force acting over an infinitesimally 
small area at a position £, referred to as a “source point.” A mathematical function 
used to describe the point force is the Dirac delta function described by Paris [6]. 
Since the problems of interest in this study have anisotropic material properties, the 
fundamental solutions for an anisotropic medium are used to develop the required 
boundary element equations. The fundamental solutions for an anisotropic body 
were determined by Cruse using a complex Airy stress function and are reproduced 
below [4]: 
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(i,j, k = 1,2) 


and the A,, values are obtained from the solution of a 4x4 matrix of equations that 
are not presented here because of their excessive length, but may be found in [4]. 

The choice for the location of the evaluation point q given in eq. 10 is 
arbitrary; however, it is necessary to evaluate all the terms along the boundary T to 
form the BIE. A consequence of evaluating the arbitrary point c, along the 
boundary is that singularities will occur at certain points within the fundamental 
solutions given by eqs. 11 and 12. These special cases of singular boundary points 
are referred to as coincident source and field points. As the distance between the 
source and field points shrinks to zero, the singularity occurs during the numerical 
integration of these boundary integrals. Therefore, the following equation results 
from the evaluation of eq. 10 along the boundary and is the integral form of the BIE 



C ij (^)U/^) + | r Vi] (q, x)u j(x) • dr = J r u*. (£, x)py (x) • dr 


(13) 


where Cy constitutes a 2x2-coefficient matrix for u/ that occurs because of a 
singularity on the boundary and i and j can have the values of one or two. The 
precise values of cij are determined by evaluating the boundary integrals in eq. 13 
by using Cauchy principal values that are dependent on the local boundary 
geometry in the neighborhood of the point under evaluation. However, explicit 
evaluation of the terms in the Cy matrix is not necessary since it was shown by 
Brebbia that the effects of a unit rigid-body displacement might indirectly calculate 
them [7]. Equation 13 may also be written using matrix notation as, 

Cu + j^ P*u-dr=J U*p dT (14) 

Using eq. 14 and the boundary conditions from eqs. 6 and 7, the unknown 
displacements and tractions on the boundary are determined. 


Interior Displacements, Strains, and Stresses 


The interior (boundary excluded) displacements, strains, and stresses are all 
determined by using Somigliana’s identity once the unknown boundary values have 
been determined. The interior displacements are calculated directly from eq. 10, 
while the strains and stresses require additional attention. Differentiating eq. 10, the 
strain tensor is obtained using the kinematic relation of eq. 2 
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where the derivatives of the fundamental solutions in eq. (15) have been evaluated 
by Cruse [4] and are reproduced below as 


S ju = -2 • 91e \ R„Q, ;i (iqfq - n 2 ) + R /2 Q, :2 (M^rq - n 2 ) 
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where the indices i, j, and l can have values equal to one or two, Py and Qij are as 
defined in eqs. 11 and 12, and Rj/e = 1 and R 2k = \Xk for k = 1,2. Thus, the interior 
strains are obtained directly from eqs. 16-18. The interior stresses are then 
determined by substituting the calculated strains into the constitutive relationship 
given by eq. 4. 

Solution Procedure 

The standard BEM solves complex problems by discretizing the boundary into 
elements and nodes, and this approach allows piecewise solution of a BIE like that 
given in eq. 14. To allow a better numerical treatment of a BIE that accurately 
represents curved boundaries, eq. 14 is rewritten in a discretized element format as 

N N 

°k u k (i9) 

k=\ k=\ 


where uk and p/e are the displacements and tractions, respectively, for node “k” and 
N is the total number of nodes for the problem. The variable l represents the 
arbitrary source point and the H and G matrices are written as 


and 


NEL 

XL p> dr 

L=1 L 

(20) 

NEL 

Ej r “> dr 

L=1 L 

(21) 


where NEL represents the total number of elements comprising the discretized 
region, T L is the boundary curvature of the L th element, and (f> is an array containing 
the assumed shape functions. Since uk in eq. 19 contain unknown displacements, 
eq. 19 is simplified by combining the ck and H Ik matrices for the two unknown 
displacement vectors in the following manner, 

when l = k, n lk =n lk +c k and 
when l*k, H lk =H lk . 


Simplifying eq. 19, 

N N 

( 22 ) 

k=l k-\ 

or in matrix form 


Hu = Gp 


(23) 



which represents all the contributions by the source point l into a global system of 
equations. Equation 23 will produce a system of equations in which known and 
unknown values are scattered between the two vectors, u and p. The known and 
unknown values are redistributed in a manner that places all known values in one 
vector and all unknown values in the other vector, which yields a desirable system 
of equations in the standard form Ax=b. 

Singular Integrals 

The method presented in this paper involves numerical integration of the boundary 
integral equation as opposed to exact solution of the individual terms in eqs. 20 and 
21. All the integrals containing non-singular functions in this paper are evaluated 
using the 8-point, regular form of Gaussian quadrature given by 

,^ 4 ^ (24) 

i=l 

where Wi are the weighting functions and q, are the integration points. In the cases 
where the integrand has a logarithmic singularity, it is necessary to manipulate the 
integrand to obtain a form that may be evaluated using the logarithmic Gaussian 
quadrature or a combination of the regular and the logarithmic Gaussian quadrature 
forms. The logarithmic Gaussian quadrature form is written as, 

I= L Zn (i) Ax)dxs i w ^) 

where the weighting function Wi and the integration point q, are not the same as 
those for the regular Gaussian quadrature. 

An important aspect of the BEM development is the relationship between an 
arbitrary boundary element and a general source point. The orientation of the 
boundary element, its distance from the source point, and the locations of the nodal 
positions are illustrated in Figure 2a. The element represented in Figure 2 is 



Figure 2: Actual and mapped coordinate systems for an element; 
a) original (actual) element and b) mapped element. 



quadratic and is used to describe both the geometry and displacement values for this 
method. The element shown in Figure 2b represents the mapping of the originally 
curved, quadratic element into a transformed coordinate system. Initially, a 
description of the x and y coordinate of any point along an element is given by its 
nodal coordinates in conjunction with the corresponding shape functions, i.e., 


x = tyl x k + $2 x k+l + fa x k+2 (26) 

y = ^y *+4>2y*+i- H l>3y*+2 ( 27 ) 

where 

fc=-±5(l-5) (28) 

<b=l-5 2 (29) 

9 3 =±$(1+5) (30) 


The distance from the source point p to a field point g on the boundary is the square 
root of the square of the x and y distances between the points. However, for an 
anisotropic material, the incorporation of complex variables into the solution 
requires that the description of the distance to a field point on the boundary take the 
following form 


z k =(x g -Xp) + m.(y g -y p ) (31) 

where \\. k are the roots of the characteristic equation defined by eq. 8. 

The singularities of the BIE are determined by finding the locations where 
the source and field points are coincident in the H and G matrices defined by eq. 
23. Cursory examination of the H and G matrices reveals that the contributions 
from the source point, i.e. the singular terms, only occur along the diagonal of each 
matrix. The location of these singular terms along the diagonals is an important 
characteristic that will allow the singularities of the H matrix in eq. 20 to be dealt 
with directly. First, recall from the section on formulating the BIE that the diagonal 
of the H matrix was determined by combining the a and H„; terms because of the 
effects of a unit rigid-body displacement on the domain Q.. In particular, the right 
hand side of eq. 22 must equal zero for a unit rigid-body displacement, which 
allows the diagonal elements of the H matrix to be expressed in terms of the off- 
diagonal elements of the H matrix. Therefore, the explicit evaluation of all the 
singular terms in the H matrix has been avoided by assuming a unit rigid-body 
displacement. The G matrix, in contrast, requires direct evaluation of each singular 
integral along its diagonal. The elements of the G matrix are comprised of the 
shape-function matrix (f> and the fundamental solution u . The shape-function 
matrix is simply a function of the Gaussian coordinate 5 an d, consequently, 
contains no singularities. Therefore, the only function that contains a singularity is 
the fundamental solution for the displacements given by eq. 11. In order to deal 
with the singularities in the G matrix, the functions must be manipulated to allow 
Gaussian integration of each element. In addition, three special cases of coincident 
source and field points exist which produce a singularity. An example for each of 




Figure 3: A general boundary element and the three special cases of coincident source and field 

points; a) l=k, b) l=k+ 1, and c) l=k+ 2. 


these cases is illustrated in Figure 3 where the source point coincides with some 
field point k; such that l=k, l=k+ 1, or l=k+ 2. 

The singular integrands of the G matrix must be rewritten to deal with the 
singularity at the coincident source and field points. The values of the natural 
logarithm function of the complex variables Zi and Z 2 in eq. 11 are the only terms 
that are not a constant, and thus contain the terms that become singular. The first 
step towards managing the singular functions of the G matrix involves rewriting the 
expression for the distance between the source and the field points. By expanding 
the complex variable zk, a more convenient form of the variable is obtained that 
both fits Gaussian quadrature and allows isolation of the singularity. This change 
of format is shown in eq. 31 as the distance from the source point to the field point 
in terms of rectangular coordinates is more conveniently written using polar 
coordinates as 


ln(z k )=ln(v k e lSk ) 


(32) 


where 
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Now, consider the case when the source and the field points are coincident ( l=k , in 
Figure 3a). The singular integral is of the form 


Jzrc(re ie )<|> r ds j= 1,2,3. (33) 

0 

Separating the kernel of eq. 33, rewriting the integrals in Gaussian integration 
format, and setting them equal to an arbitrary value Ii for convenience, such that 

+1 +1 

I| ={ ln[v(mi © det[J] • + J ie ©<|> 2 © det[J] • d^ 

-1 -1 


( 34 ) 



which can be rewritten as Ii = I2 + I3. The second integral, I3, is nonsingular since it 
approaches zero as the radius goes to zero, so I3 in equation 34 can be evaluated by 
simply using regular Gaussian quadrature. In order to evaluate the first integral, I2, 
a change of limits is required because of the singularity. First, the limits of P are 
changed to 0 to 1 (from -1 to +1) 

1 

I 2 = J/re[rCn)M» 2 (Ti)det[J]-2-dTi ( 35 ) 
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and then the numerator and denominator of the kernel are multiplied by p and 
separated as shown below 
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which can again be rewritten as I2 = I4 + I5. Thus, the problem of the singularity has 
been resolved in F by splitting the integral into two integrals that are more 
manageable. The first integral I4 can be shown to be non-singular using L'Hopital's 
rule which enables it to be integrated using the regular Gaussian quadrature. The 
second integral, I5, has been conveniently manipulated into the logarithmic 
Gaussian quadrature format, such that 
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( 37 ) 


Returning to I4, the limits of integration need to be adjusted to facilitate the use of 
regular Gaussian quadrature, so 
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Finally, for the case of coincident source and field point, when l=k , we can write the 
integral I] as, 
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where N rg are evaluation points used for the regular form of Gaussian quadrature 
and Ni g are evaluation points used for the logarithmic form of Gaussian quadrature. 
Equation 39 does not contain any singularities and may be evaluated at any point 
under consideration. Also, note that a sequence of coordinate mappings was used to 
determine eq. 39. Therefore, it is very important to remember the critical 
information; the radius r(^), the shape function (j y, and the Jacobian J must always 
be calculated in terms of the original coordinate The process of removing 
singularities for the other two cases of coincident source and field points, when 
l=k + 1 and l=k+ 2, are very similar to the previous derivation and the complete 
derivations are available in earlier work by Smeltzer [10]. 


NUMERICAL RESULTS 

The BEM presented in the previous section was evaluated by using a quasi- 
isotropic plate with an elliptical cutout. The quasi-isotropic plate configuration 
discussed in the present paper was chosen for convenience, and essentially 
demonstrates the full capabilities of the curved, quadratic boundary element. The 
plate under consideration is square with an edge length of 10 inches and was loaded 
by subjecting two of the straight edges of the plate to uniform tension in the 
direction of the xi-axis. The elliptical cutout is at the center of the plate and has a 
major axis of 0.1 inches and a minor axis of 0.05 inches where the minor axis is 
aligned with the x r axis. A half-plate model of the plate with the appropriate 
symmetry conditions was used with 54 quadratic boundary elements, and only eight 
of those elements were used to model the boundary of the elliptical cutout. The 
effective, laminate material properties for the [0,±45] s plate are given in Table 1. In 
addition, the subscripts 1 and 2 for the material properties given in Table 1 refer to 
the coordinate directions of the structure and not the principal material coordinate 
directions of a laminate ply. 

The analytical solution for an infinite plate is given by Lekhnitskii and is 
plotted in Figure 4 along with the BEM solution [8]. Correlation between the 
analytical solution and the corresponding BEM solution is excellent, with an error 
of less than 1%. In addition, the figure shows a stress concentration at the edge of 
the elliptical cutout equal to 3.74 times the far-field traction P. 

TABLE 1: BASIC MATERIAL PROPERTIES FOR THE COMPOSITE PLATE 


Laminate 

Roots of Characteristic 
Equation (8) 

Effective, Laminate 
Material Properties 

[0,±45] s 

lii = 0.7921 +0.9659i 
p 2 = -0.7921 +0.9659i 

Ej = 7.706 x 10 6 psi 
E 2 = 3.164 x 10 6 psi 
v 12 = 0.794 
Gi 2 = 3.504 x 10 6 psi 





Figure 4: The stress distribution along the positive x 2 -axis for a [0,±45] s 
composite laminate with an elliptical cutout. 


CONCLUDING REMARKS 

A method has been presented in the present paper that details the development of a 
novel, curved, quadratic boundary element for analyzing structures composed of 
anisotropic materials in a state of plane stress or plane strain. Although the standard 
boundary element method (BEM) equations for elastostatics were employed in the 
development of the curved boundary element, elaborate singular integrals were 
derived because of the nature of the complex fundamental solutions for an 
anisotropic material. The treatment of these singular integrals was a key component 
of the boundary element derivation presented herein, and was complicated further 
by allowing the boundary element to be curved. The complexity of the singular 
integrals was handled by performing a sequence of coordinate mappings on the 
curved elements that resulted in an expression that was amenable to Gaussian 
quadrature. Thus, the curved elements were mapped to a new coordinate system 
that facilitated removal of the singularities and allowed for accurate and efficient 




evaluation of the integrals. The basic framework and solution technique for the 
curved BEM was illustrated and a comparison of the results for a curved, quadratic 
boundary element with those for a classic problem in engineering mechanics 
showed differences of less than 1%. 
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